################################################################################
########################     Analysis and Plots     ############################
################################################################################

library(Matching)
library(tidyverse)

################################################################################

# Input is virginclip_pop.RData
# I omit the command to read in the dataset as it will need a directory.
# save.image("~/Dropbox/PHD/Indonesia_GPS/results/dryland-1618.RData")
rm(list=ls())
load("~/data/virginclip_pop.RData")

################################################################################

# 1. No trimming on forest cover in 2005

virginclip_pop <- virginclip_pop %>%
  filter(peat_depth_GW==0)

# outcome vars

virginclip_pop$did2016 <- virginclip_pop$forestha_2016 - virginclip_pop$forestha_2011
#virginclip_pop$did2005 <- virginclip_pop$forestha_2010 - virginclip_pop$forestha_2005

virginclip_pop$did2017 <- virginclip_pop$forestha_2017 - virginclip_pop$forestha_2011
#virginclip_pop$did2004 <- virginclip_pop$forestha_2010 - virginclip_pop$forestha_2004

virginclip_pop$did2018 <- virginclip_pop$forestha_2018 - virginclip_pop$forestha_2011
#virginclip_pop$did2003 <- virginclip_pop$forestha_2010 - virginclip_pop$forestha_2003

t1 = 2010
t11 = 2011
t0 = 2004

t2 = 2016
t3 = 2017
t4 = 2018

virginclip_pop$didpre = virginclip_pop$forestha_2010 - virginclip_pop$forestha_2004

virginclip_pop$ddd2016 <- virginclip_pop$did2016 - ((t2-t11)/(t1-t0))*virginclip_pop$didpre
virginclip_pop$ddd2017 <- virginclip_pop$did2017 - ((t3-t11)/(t1-t0))*virginclip_pop$didpre
virginclip_pop$ddd2018 <- virginclip_pop$did2018 - ((t4-t11)/(t1-t0))*virginclip_pop$didpre

# ps model
ps_virginclip_pop <- glm(moratorium ~ forestha_2005 + forestha_2006 + forestha_2007 + 
                       forestha_2008 + forestha_2009 + forestha_2010 +  elev + slope + 
                       distroads + dist_cities + dist_palms_km + dist_timber_km + dist_logging_km +
                       agc_bcc_mean + pop_2005 + pop_2010, family = binomial, data = virginclip_pop)
summary(ps_virginclip_pop)

virginclip_pop$pscore<- ps_virginclip_pop$fitted

###### DID REGRESSIONS ######

PSM_DID_virginclip_pop2016 <- Match( Y = virginclip_pop$did2016, Tr = virginclip_pop$moratorium,
                                 X = virginclip_pop$pscore, estimand = "ATT", M = 1, ties = T,
                                 replace = T, distance.tolerance = 1e-15, caliper = 0.001)

summary(PSM_DID_virginclip_pop2016)




PSM_DID_virginclip_pop2017 <- Match( Y = virginclip_pop$did2017, Tr = virginclip_pop$moratorium,
                                 X = virginclip_pop$pscore, estimand = "ATT", M = 1, ties = T,
                                 replace = T, distance.tolerance = 1e-15, caliper = 0.001)

summary(PSM_DID_virginclip_pop2017)

PSM_DID_virginclip_pop2018 <- Match( Y = virginclip_pop$did2018, Tr = virginclip_pop$moratorium,
                                 X = virginclip_pop$pscore, estimand = "ATT", M = 1, ties = T,
                                 replace = T, distance.tolerance = 1e-15, caliper = 0.001)

summary(PSM_DID_virginclip_pop2018)



###### DDD REGRESSIONS ######


PSM_DDD_virginclip_pop2016 <- Match( Y = virginclip_pop$ddd2016, Tr = virginclip_pop$moratorium,
                                 X = virginclip_pop$pscore, estimand = "ATT", M = 1, ties = T,
                                 replace = T, distance.tolerance = 1e-15, caliper = 0.001)

summary(PSM_DDD_virginclip_pop2016)

PSM_DDD_virginclip_pop2017 <- Match( Y = virginclip_pop$ddd2017, Tr = virginclip_pop$moratorium,
                                 X = virginclip_pop$pscore, estimand = "ATT", M = 1, ties = T,
                                 replace = T, distance.tolerance = 1e-15, caliper = 0.001)

summary(PSM_DDD_virginclip_pop2017)


PSM_DDD_virginclip_pop2018 <- Match( Y = virginclip_pop$ddd2018, Tr = virginclip_pop$moratorium,
                                 X = virginclip_pop$pscore, estimand = "ATT", M = 1, ties = T,
                                 replace = T, distance.tolerance = 1e-15, caliper = 0.001)

summary(PSM_DDD_virginclip_pop2018)


rm(virginclip_pop, ps_virginclip_pop)

save.image("~/results/dryland-cf-pop-161718-0001.RData")


### peatland ###

rm(list=ls())
load("~/data/virginclip_pop.RData")

################################################################################

# 1. No trimming on forest cover in 2005

virginclip_pop <- virginclip_pop %>%
  filter(peat_depth_GW>0)

# outcome vars

virginclip_pop$did2016 <- virginclip_pop$forestha_2016 - virginclip_pop$forestha_2011
#virginclip_pop$did2005 <- virginclip_pop$forestha_2010 - virginclip_pop$forestha_2005

virginclip_pop$did2017 <- virginclip_pop$forestha_2017 - virginclip_pop$forestha_2011
#virginclip_pop$did2004 <- virginclip_pop$forestha_2010 - virginclip_pop$forestha_2004


t1 = 2010
t11 = 2011
t0 = 2004

t2 = 2016
t3 = 2017

virginclip_pop$didpre = virginclip_pop$forestha_2010 - virginclip_pop$forestha_2004

virginclip_pop$ddd2016 <- virginclip_pop$did2016 - ((t2-t11)/(t1-t0))*virginclip_pop$didpre
virginclip_pop$ddd2017 <- virginclip_pop$did2017 - ((t3-t11)/(t1-t0))*virginclip_pop$didpre


# ps model
ps_virginclip_pop <- glm(moratorium ~ forestha_2005 + forestha_2006 + forestha_2007 + 
                       forestha_2008 + forestha_2009 + forestha_2010 +  elev + slope + 
                       distroads + dist_cities + dist_palms_km + dist_timber_km + dist_logging_km +
                       agc_bcc_mean + pop_2005 + pop_2010, family = binomial, data = virginclip_pop)
summary(ps_virginclip_pop)

virginclip_pop$pscore<- ps_virginclip_pop$fitted

# DDD models


PSM_DDD_virginclip_pop2016 <- Match( Y = virginclip_pop$ddd2016, Tr = virginclip_pop$moratorium,
                                 X = virginclip_pop$pscore, estimand = "ATT", M = 1, ties = T,
                                 replace = T, distance.tolerance = 1e-15, caliper = 0.001)

summary(PSM_DDD_virginclip_pop2016)

PSM_DDD_virginclip_pop2017 <- Match( Y = virginclip_pop$ddd2017, Tr = virginclip_pop$moratorium,
                                 X = virginclip_pop$pscore, estimand = "ATT", M = 1, ties = T,
                                 replace = T, distance.tolerance = 1e-15, caliper = 0.001)

summary(PSM_DDD_virginclip_pop2017)



# remove data and save results
rm(virginclip_pop, ps_virginclip_pop)

save.image("~/results/peatland-cf-pop-1617-0001.RData")